#load required packages
library(EasyABC)
library(abcrf)
library(ggplot2)

#load required data
load("sum_stats.RData")

#set up function for ABC with bias as a parameter, using a modified version of the transmission function in HERAChp.KandlerCrema that uses less memory
transmission_model <- function(parameter){
  source("hpcc/modified_herachp_more_output.R")
  set.seed(parameter[1])
  
  samples <- transmission(N = 729, timesteps = 232, mu = 0.03672527, warmUp = 200, top = 142, bias = parameter[2])
  return(c(samples$x, mean(samples$obs.div), samples$obs.div, colMeans(samples$zMatrix)))
  
  rm(samples)
}

#simulate conformity model
priors <- list(c("unif", -0.2, 0))
abc_conformity <- ABC_rejection(model = transmission_model, prior = priors, nb_simul = 50000, use_seed = TRUE, n_cluster = 11)

#simulate novelty model
priors <- list(c("unif", 0, 0.2))
abc_novelty <- ABC_rejection(model = transmission_model, prior = priors, nb_simul = 50000, use_seed = TRUE, n_cluster = 11)

#set up function for ABC without bias as a parameter, using a modified version of the transmission function in HERAChp.KandlerCrema that uses less memory
transmission_model <- function(parameter){
  source("hpcc/modified_herachp_more_output.R")
  set.seed(parameter[1])
  
  samples <- transmission(N = 729, timesteps = 232, mu = 0.03672527, warmUp = 200, top = 142, bias = 0)
  return(c(samples$x, mean(samples$obs.div), samples$obs.div, colMeans(samples$zMatrix)))
  
  rm(samples)
}

#simulate neutrality model
abc_neutrality <- matrix(data = NA, nrow = 50000, ncol = 176)
i <- 1
while(i <= 50000){
  abc_neutrality[i, ] <- transmission_model(i)
  i <- i + 1
}

#generate data frame for random forest
sumsta <- rbind(abc_conformity$stats, abc_novelty$stats, abc_neutrality)
sumsta <- as.data.frame(sumsta)
modindex <- as.vector(c(rep("1", 50000), rep("2", 50000), rep("3", 50000)))
abcrf_data <- data.frame(modindex, sumsta)

#run random forest
abcrf <- abcrf(formula = modindex ~ ., data = abcrf_data, ntree = 1000, sampsize = 150000, paral = TRUE, ncores = 11)

sum_stats_obs <- data.frame(t(sum_stats))
colnames(sum_stats_obs) <- colnames(abcrf_data)[-1]

#tally votes
abcrf_predict <- predict(abcrf, sum_stats_obs, abcrf_data, paral = TRUE, ncores = 11)
